
#include <msGeneralizedCoordinates.h>
#include <msParser.h>
#include <msZMat.h>
#include <msSystem.h>
#include <msKineticOperator.h>
#include <msLagrangian.h>
#include <muParser.h>
#include <msQuantumMotion1d.h>
#include <msSolverLagrangian.h>


int main(int argc, char **argv) {
 
 using namespace atomism;

 // Test object creation
 boost::shared_ptr<msUnitsManager> unitMng = 	msUnitsManager::New();

 unitMng->set( "Angstrom Degree amu kcal/mol ps" );

 boost::shared_ptr<msScalarVariable> q0 = 	msScalarVariable::New("Degree");
 q0->set(0,0,360,1,1,"Degree");
 q0->setId("q0");


 boost::shared_ptr<msGeneralizedCoordinates> gcoors = 
						msGeneralizedCoordinates::New(unitMng);
 gcoors->addVar(q0);

 boost::shared_ptr<msParser> kinfct14 = 	msParser::New("Angstrom");
 kinfct14->setCoordinates(gcoors);
 kinfct14->setExpression("q0");
 kinfct14->setId("f=q0");

 boost::shared_ptr<msParser> kinfct17 = 	msParser::New("Angstrom");
 kinfct17->setCoordinates(gcoors);
 kinfct17->setExpression("q0+180");
 kinfct17->setId("f=q0+180");

 boost::shared_ptr<msZMat> zmt =  		msZMat::New(unitMng);
 zmt->set("C C 0 1.34 H 0 1.09 1 121.8 H 0 1.09 2 116.6 1 180 C 1 1.46 0 124.3 2 180 H 1 1.09 0 119.4 2 0 C 4 1.34 1 124.3 0 -180 H 4 1.09 1 116.2 0 0 H 6 1.09 4 121.5 1 0 H 6 1.09 4 121.8 1 180");
 zmt->setId("zmat");

 
 boost::shared_ptr<msSystem> system =  		msSystem::New(unitMng);
 system->addEntity(zmt);
 
 boost::shared_ptr<msKineticOperator> kinop =  	msKineticOperator::New(unitMng);
 kinop->set(gcoors,system);
 kinop->addDynamicDof("zmat",14,kinfct14);
 kinop->addDynamicDof("zmat",17,kinfct17);
 
 boost::shared_ptr<msParser> epot = msParser::New("kcal/mol");
 epot->setCoordinates(gcoors);
 epot->setExpression("5*(1+cos(q0/180*3.14))");
 epot->setId("f=epot(q0)");

 boost::shared_ptr<msLagrangian> motion = 	msLagrangian::New(unitMng);
 motion->setEqOfMotion(kinop,epot);

 boost::shared_ptr<msQuantumMotion1d> qmotion = msQuantumMotion1d::New(unitMng);
 qmotion->setEqOfMotion(kinop,epot);
 boost::shared_ptr<msVectorFit1d> Kineff = qmotion->computeEffectiveKineticFunction();
 qmotion->computeEigenValues();
 boost::shared_ptr<msVectorFit1d> ES0    = qmotion->eigenvector(0); 
 try{
     qmotion->getParameters()->setParameter("Emax",	20,		*unitMng);
     qmotion->getParameters()->setParameter("BasisSize",201,		*unitMng);
     qmotion->getParameters()->setParameter("A joke","raise an error",	*unitMng);
     }
 catch(msError& e)
     {
      e.addContext("this was just a test for parameter exception");
      cout<<e.what();
      fprintf(stdout,"exception caught for bad parameters.");
      }

 boost::shared_ptr<atomism::msSolverLagrangian> solver = atomism::msSolverLagrangian::New();
 //----------------------------------------------------------

 // test for segmetation fault in print functions 
 unitMng->printUnityMap();

 cout<<*unitMng<<endl;
 cout<<*q0<<endl;
 cout<<*gcoors<<endl;
 cout<<*kinfct14<<endl;
 cout<<*zmt<<endl;
 cout<<*system<<endl;
 cout<<*kinop<<endl;
 cout<<*motion<<endl;
 
 cout<<*qmotion<<endl;
 cout<<*zmt<<endl;
 cout<<*Kineff<<endl;
 cout<<*ES0<<endl;
 cout<<*solver<<endl;
 motion->getParameters()->listWhoUsingMe();
 unitMng->listWhoUsingMe();
 motion->getSampler()->print(cout);
 

 // test calculs 

 if( ES0->max() != 0.0741 ) {}//{ throw; };

 solver->dynamic(*motion,10);
 fprintf(stdout,"Efinal=%g ; Tfinal=%g ; Lfinal=%g ; ",	motion->E(),motion->T(),motion->L());
 fprintf(stdout,"kin=%g ; ",				kinop->kineticFunction());
 msKineticWork work = kinop->move(0);
 fprintf(stdout,"move: dJ=%g , dP=%g ; ",		work.dJ.norm(), work.dP.norm());
 fprintf(stdout,"eff. kin(300)=%g ; ",			qmotion->effectiveKineticFunction(300,*unitMng)); 
 fprintf(stdout,"efF. kin(300) from eff. kin. fct =%g ; ",(*Kineff)(300) );
 fprintf(stdout,"Qclass(300)=%g ; ",			motion->Q(300));
 fprintf(stdout,"Qquant(300)=%g ; ",			qmotion->Q(300));
 fprintf(stdout,"QquantFromDos(300)=%g ; ",		qmotion->QfromDOS(300));
 fprintf(stdout,"DOSclass(300)=%g ; ",			motion->DOS(5));
 fprintf(stdout,"DOSquant(300)=%g\n",			motion->DOS(5));
 fprintf(stdout,"SOSquant(300)=%g\n",			motion->SOS(5));
 fprintf(stdout,"Qquant(300,kinEffFct)=%g ; ",		qmotion->Q(300,*Kineff));  

 
 // test the tree mapper : add child, remove child etc

 qmotion->addChild( ((msMotion*)qmotion.get())->DOS() , "dos computed" );

 boost::shared_ptr<msUnitsManager>    unitMng2 	= boost::static_pointer_cast<msUnitsManager>(    unitMng->deepCopy() );
 boost::shared_ptr<msScalarVariable>    q2 	= boost::static_pointer_cast<msScalarVariable>(  q0->deepCopy()      );
 boost::shared_ptr<msGeneralizedCoordinates>    gcoors2 = boost::static_pointer_cast<msGeneralizedCoordinates>(  gcoors->deepCopy()      );
 boost::shared_ptr<msParser>    epot2 		= boost::static_pointer_cast<msParser>(  epot->deepCopy()      );
 boost::shared_ptr<msZMat>    zmt2 		= boost::static_pointer_cast<msZMat>(  zmt->deepCopy()      );
 boost::shared_ptr<msSystem>    system2 	= boost::static_pointer_cast<msSystem>(  system->deepCopy()      );
 boost::shared_ptr<msKineticOperator>    kinop2 = boost::static_pointer_cast<msKineticOperator>(  kinop->deepCopy()      );
 boost::shared_ptr<msLagrangian> motion2 	= boost::static_pointer_cast<msLagrangian>( motion->deepCopy() );
 boost::shared_ptr<msQuantumMotion1d> qmotion2 	= boost::static_pointer_cast<msQuantumMotion1d>( qmotion->deepCopy() );


// std::vector<double> a;
 //std::cout<<a<<std::endl;
/* cout<<"SOSclass(5)="<<motion->SOS(5)<<endl;
 cout<<"SOSquant(5)="<<qmotion->SOS(5)<<endl;*/

 return 0;
 }


